Studies on optimizing potential energy functions for maximal intrinsic 

hyperpolarizability 
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We use numerical optimization to study the properties of (1) the class of one-dimensional po- 
tential energy functions and (2) systems of point charges in two-dimensions that yield the largest 
hyperpolarizabilities, which we find to be within 30% of the fundamental limit. We investigate the 
character of the potential energy functions and resulting wavefunctions and find that a broad range 
of potentials yield the same intrinsic hyperpolarizability ceiling of 0.709. 



I. INTRODUCTION 

Materials with large nonlinear-optical suscep- 
tibilities are central for optical applications such 
as telecommunications, [l| three-dimensional nano- 
photolithography, 0, Q and making new materials Q for 
novel cancer therapies. || The fact that quantum calcu- 
lations show that there is a limit to the nonlinear-optical 
response d, 0, H, S OH EH is both interesting from the 
basic science perspective; and, provides a target for 
making optimized materials. In this work, we focus 
on the second-order susceptibility and the underlying 
molecular hyperpolarizability, which is the basis of 
electro-optic switches and frequency doublers. 

The fundamental limit of the off-resonance hyperpo- 
larizability is given by,Q 
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where N is the number of electrons and Eiq the energy 
difference between the first excited state and the ground 
state, Eiq = E\ — Eq. Using Equation [TJ we can de- 
fine the off-resonant intrinsic hyperpolarizability, /3j„t , as 
the ratio of the actual hyperpolarizability (measured or 
calculated), (3, to the fundamental limit, 
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We note that since the dispersion of the fundamental 
limit of (3 is also known, [12| it is possible to calculate 
the intrinsic hyperpolarizability at any wavelength. In 
the present work, we treat only the zero- frequency limit. 

Until recently, the largest nonlinear susceptibilities of 
the best molecules fell short of the fundamental limit by a 
factor of 10 3 / 2 , 03, EH so the very best molecules had 
a value of Pint = 0.03. Since a Sum-Over-States (SOS) 
calculation of the hyperpolarizability 15] using the ana- 
lytical wavefunctions of the clipped harmonic oscillator 
yields a value /3j„t = 0.57, the factor-of-thirty gap is 
not of a fundamental nature. Indeed, recently, it was re- 
ported that a new molecule with asymmetric conjugation 
modulation has a measured value of /3j n t = 0.048. [lq 

To investigate how one might make molecules with 
a larger intrinsic hyperpolarizability, Zhou and cowork- 
ers used a numerical optimization process where a trial 



potential energy function is entered as an input, and 
the code iteratively deforms the potential energy func- 
tion until the intrinsic hyperpolarizability, calculated 
from the resulting wavefunctions, converges to a local 
maximum. [l7| In this work, a hyperbolic tangent func- 
tion was used as the starting potential due to the fact that 
it is both asymmetric yet relatively flat away from the ori- 
gin. This calculation was one-dimensional and included 
only one electron, so electron correlation effects were 
ignored. Furthermore, the intrinsic hyperpolarizability 
was calculated using the new dipole-free sum-over-states 
expression [l8| and only 15 excited states were included. 
The resulting optimized potential energy function showed 
strong oscillations, which served to separate the spatial 
overlap between the energy eigenfunctions. This led Zhou 
and coworkers to propose that modulated conjugation 
in the bridge between donor and acceptor ends of such 
molecules may be a new paradigm for making molecules 
with higher intrinsic hyperpolarizability. fl7| 

Based on this paradigm, Perez Moreno reported mea- 
surements of a class of chromophobes with varying de- 
gree of modulated conjugation. [16( The best measured 
intrinsic hyperpolarizability was (3i n t — 0.048, about 50% 
larger than the best previously-reported. Given the mod- 
est degree of conjugation modulation for this molecule, 
this new paradigm shows promise for further improve- 
ments. 

In the present work, we extend Zhou's calculations to 
a larger set of starting potentials. To circumvent trun- 
cation problems associated with sum-over-states calcula- 
tions, we instead determine the hyperpolarizability using 
a finite difference technique. The optimization procedure 
is then applied to this non-perturbative hyperpolarizabil- 
ity. 

To study the effects of geometry on the hyperpo- 
larizability, Kuzyk and Watkins calculated the hyper- 
polarizability of various arrangements of point charges, 
representing nuclei, in two-dimensions using a two- 
dimensional Coulomb potential. In the present con- 
tribution, we apply our numerical optimization technique 
to determine the arrangement and charges of the nuclei 
in a planar molecule that maximize the intrinsic hyper- 
polarizability. 
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II. THEORY 

In our previous work, we used a finite-state SOS model 
of the hyperpolarizability that derives from perturbation 
theory (we used both the standard Orr and Ward SOS ex- 
pression, /SsoSiEU an d the newer dipole free expression, 
/3_d_f[3). The use of a finite number of states in lieu of 
the full infinite sums can result in inaccuracies, so, in the 
present work, we use the non-perturbative approach, as 
follows. We begin by solving the 1-d Schrodinger Equa- 
tion on the interval a < x < b for the ground state wave- 
function ip(x, E) of an electron in a potential well defined 
by V(x) and in the presence of a static electric field, E, 
that adds to the potential SV = —exE. From this, the 
off-resonant hyperpolarizability is calculated with numer- 
ical differentiation, i.e. using finite differences, yielding 
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Equation [3] is evaluated using the standard second-order 
approximation to the second derivative: 
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with several h values ho, ho/ 5, ho/ 25, ... . We then refine 
these values by Richardson extrapolation (2(i| and obtain 
our estimate from the two closest extrapolated values. 

Our computational mesh consists of 200 quadratic fi- 
nite elements with a total of 399 degrees of freedom. The 
potential energy function is a cubic spline with 40 degrees 
of freedom. Thus the numerical calculations in regions 
where the potential function is represented by 3 points 
in the spline are covered by 15 elements with a total of 
about 30 degrees of freedom. 

Calculating p int from Equations[31[2]and[T]for a specific 
potential, we use an optimization algorithm that contin- 
uously varies the potential in a way that maximizes Pint- 
We also compute the matrix [l7l l2q ] 



and energies are exact. We have found that the values 
of Tmp are small for exact wavefunctions when m < N/2 
and p < N/2. So, when evaluating the r matrix to 
test our calculations, we consider only the components 
T W 

' m<N/2,p<N/2- 

We observe that when using more than about 40 states 
in SOS calculations of the hyperpolarizability only a 
marginal increase of accuracy results when the poten- 
tial energy function is parameterized with 400 degrees of 
freedom. So, to ensure overkill, we use 80 states when 
calculating the r matrix or the hyperpolarizability with 
an SOS expression so that truncation errors are kept to 
a minimum. Since the hyperpolarizability depends crit- 
ically on the transition dipole moment from the ground 
state to the excited states, we use the value of Tqq as 
one important test of the accuracy of the calculated wave- 
functions. Additionally, we use the standard deviation of 
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which quantifies, on average, how well the sum rules are 
obeyed in aggregate, making a broader test of the 

accuracy of a large set of wavefunctions. 

Our code is written in MATLAB. For each trial po- 
tential we use a quadratic finite element method [2l[ to 
approximate the Schrodinger eigenvalue problem and the 
implicitly restarted Arnoldi method [22| to compute the 
wave functions and energy levels. T o op timize (3 we use 
the Nelder-Mead simplex algorithm |23| . 

As described in our previous work, [17j we perform op- 
timization, but this time using the exact intrinsic hyper- 
polarizability (3 — Pnp/Pmax, where (3m ax is the fun- 
damental limit of the hyperpolarizability, which is pro- 
portional to E^ . To determine Eiq = E\ — Eq, we also 
calculate the first excited state energy E\. 
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where x" ax is the magnitude of the fundamental limit of 
the position matrix element xio for a one electron system, 
and is given by, 
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Each matrix element of r"', indexed by m and p, is 
a measure of how well the (m,p) sum rule is obeyed 
when truncated to N states. If the sum rules are ex- 
actly obeyed, Tmp^ — for all m and p. We note that if 
the sum rules are truncated to an N-state model, the sum 
rules indexed by a large value of m or p (i.e. m,p ~ N) 
are disobeyed even when the position matrix elements 



III. RESULTS AND DISCUSSIONS 

Figure [1] shows an example of the optimized poten- 
tial energy function after 7,000 iterations when starting 
with the potential V(x) = and optimizing the non- 
perturbative intrinsic hyperpolarizability /?jv p / (3m ax as 
calculated with Equation [31 Also shown are the eigen- 
functions of the first 15 states computed from the opti- 
mized potential. First, we note that the potential energy 
function shows the same kinds of wiggles as in our original 
paper, [17] though not of sufficient amplitude to localize 
the wavefunctions. 

For the starting potentials we have investigated, our 
results fall into two broad classes. In the first, three 
common features are: (1) The best intrinsic hyperpolar- 
izabilities are near /3j„ t = 0.71; (2) the best potentials 
have a series of wiggles; and (3) the systems behave as a 
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FIG. 1: Optimized potential energy function and first 15 
wavefunctions after 7,000 iterations. Starting potential is 
V(x) — 0, using the non-perturbative hyperpolarizability for 
optimization. 



limited-state model. In the second class of starting po- 
tentials, (2) the wiggles are much less pronounced and (3) 
more states contribute evenly. Figure [T] is an example of 
a Class II potential. However, in both classes, the max- 
imum calculated intrinsic hyperpolarizability appears to 
be bounded by 0i n t = 0.71. Using the set of potentials 
from both classes that lead to optimized (3n p / Pmax , we 
calculate the lowest 80 eigenfunctions and eigenvalues, 
from which we calculate (3df and fisos- hi most cases, 
we find that the three different formulas for (3 converge 
to the same value when only the first 20 excited states 
are used (using 80 states, the three are often the same to 
within at least 4 decimal places) and Too ~ 10~ 4 , showing 
that the ground state sum rules are well obeyed. Further- 
more, the rms deviation of the r matrix when including 
40 states leads to < 0.001. 

Figure [5] shows an example of the optimized potential 
energy function when starting with the potential V(x) — 
tanhx and optimizing the exact (non-perturbative) in- 
trinsic hyperpolarizability. Also shown are the eigenfunc- 
tions of the first 15 states computed with the optimized 
potential. First, we note that the potential energy func- 
tion shows the same kinds of wiggles as in our original 
paper; [17| and only 2 excited state wavefunctions and the 
ground state are localized in the first deep well - placing 
this system in Class I. 

The observation that such potentials lead to hyper- 
polarizabilities that are near the fundamental limit mo- 
tivated Zhou and coworkers to suggest that molecules 
with modulated conjugation may have enhanced intrin- 
sic hyperpolarizabilities. 17] A molecule with a modu- 
lated conjugations bridge between the donor and accep- 
tor end was later shown to have record-high intrinsic 
hyperpolarizability. [l6j As such, this result warrants a 
more careful analysis. 

It is worthwhile to compare our present results charac- 
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FIG. 2: Optimized potential energy function and first 15 
wavefunctions after 8,000 iterations. Starting potential is 
V(x) = tanh(x), using the non-perturbative hyperpolarizabil- 
ity for optimization. 



terized by Figure [5] with our past work, [13] particularly 
for the purpose of examining the impact of the approx- 
imations used in the previous work. [17| Figure [3] shows 
the optimized potential and wavefunctions obtained by 
Zhou and coworkers using a 15-state model and opti- 
mizing the dipole-free intrinsic hyperpolarizability. Since 
only 15 states were used, the SOS expression for /3 did not 
fully converge; making the result inaccurate as suggested 
by the fact that Psos ancl Pdf did not agree. How- 
ever, since the code focused on optimizing the dipole-free 
form of (3, and Too was small when /3j ni was optimized, 
the dipole-free expression may have converged to a rea- 
sonably accurate value while the commonly-used SOS 
expression was inaccurate. Indeed, it was found that 
[3df ~ 0.72 - in contrast to our more precise present 
calculations using the non-perturbative approach, which 
yields /3jvp < 0.71. So, the fact that our more precise 
calculations, which do not rely on a sum-over states ex- 
pression, agree so well with the 15-state model suggests 
that in both cases, the limit for a one-dimensional single 
electron molecule is just over (3 w 0.7. This brute force 
calculation serves as a numerical illustration of the obser- 
vation that the limiting value of (3 is the same for an exact 
non-perturbation calculation and for a calculation that 
truncates the SOS expression, which presumedly should 
lead to large inaccuracies. [24l |25| At minimum, this result 
supports the existence of fundamental limits of nonlinear 
susceptibilities that are in line with past calculations. 

To state Zhou's approach more precisely. [l7| the cal- 
culations optimized the very special case of the intrin- 
sic hyperpolarizability for a 15 state model for a poten- 
tial energy function that is parameterized with 20 spline 
points. As such, the potential energy function can at 
most develop about 20 wiggles. As a consequence, there 
are enough degrees of freedom in the potential energy 
function to force the 15 states to be spatially well sepa- 
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FIG. 3: Optimized potential energy function using Pdf and 
first 15 wavefunctions after 7,000 iterations. Starting po- 
tential is the tanh(a;) potential. The final potential (shown 
above) we refer to as the Zhou potential. 



FIG. 4: Optimized potential energy function and first 15 
wavefunctions after 1,000 iterations. Starting potential is 
Zhou's potential, using the non-perturbative hyperpolariz- 
ability for optimization. 



TABLE I: Evolution of Zhou's Potential. S is the hyperpo- 
larizability of the starting potential using 80 states while the 
other ones are after optimization of Pnp- 



Number of 
Iterations 


Ps 


Psos Pdf Pnp 


T (8°) 

'oo 
(xlO -5 ) 


Ar (80) 

(xHT 4 ) 




1000 


0.5612 
0.5612 


0.5612 0.5607 0.5612 
0.7087 0.6682 0.7083 


11.2 
1810 


15 
40 



rated. Interestingly, after optimization, only two excited 
states overlap with ground state, allowing only these two 
states to have nonzero transition dipole moments with 
each other and the ground state - forcing the system 
into a three-level SOS model for (3df- This behavior is 
interesting in light of the three- level ansatz, which asserts 
that only three states determine the nonlinear response 
of a system when it is near the fundamental limits. 

It is interesting to compare the exact non-perturbation 
calculation, which does not depend on the excited state 
wavefunctions (Figure and Zhou's contrived system of 
15 states (Figure [3]). Both cases have wiggles and the 
wavefunctions appear to be mostly non-overlapping. So, 
for the first 15 states, the wavefunctions appear simi- 
larly localized. The situation becomes more interesting 
when 80 states are included in calculating the hyperpo- 
larizability for Zhou's potential or when the exact non- 
perturbative approach is used. The first line in Table 
U summarizes the results with Zhou's potential and 80 
states. 

First, let's focus on the sum-over-states results. 
Clearly, when 80 states are used in the calculation, it 
is impossible for the excited state wavefunctions to not 
overlap with each other, so the three-level approximation 
to (3 breaks down. According to the three-level ansatz, 
we would expect the hyperpolarizability to get smaller. 
Indeed, the additional excited states result in a smaller 



hyperpolarizability (w 0.56). Note that the exact and 
SOS expressions agree with each other and that Tqq ^ 
and At^ 80 ) are small. 

Figure 2] shows the result after 1000 iterations, us- 
ing Zhou's potential as the starting potential and us- 
ing the non-perturbative hyperpolarizability for opti- 
mization. First, the non-perturbative hyperpolarizabil- 
ity reaches just under 0.71, but, the SOS and dipole- 
free expressions do not agree with each other. Further- 
more, both convergence metrics (r ( * 0) and Ar^ 80 )) arc 
larger than before optimization. It would appear that 
for Zhou's potential, even 80 states are not sufficient 
to characterize the nonlinear susceptibility when a sum- 
over-states expression is used (either dipole free or tradi- 
tional SOS expression - though the SOS expression agrees 
better with the non-perturbative approach). 

Interestingly, the optimized potential energy function 
still retains the wiggles and the wave functions are still 
well separated. This result is consistent with the sug- 
gestion of Zhou and coworkers that modulation of conju- 
gation may be a good design strategy for making large- 
hyperpolarizability molecules. We note that wiggles in 
the potential energy function are not required to get 
a large nonlinear-optical response; but, appears to be 
one way that Mother Nature optimizes the hyperpo- 
larizability. Since this idea has been used to identify 
molecules with experimentally measured record intrin- 
sic hyperpolarizability. [l6j the concept of modulation of 
conjugation warrants further experimental studies. 

As a case in point that non-wiggly potentials can lead 
to a large nonlinear susceptibility is the clipped harmonic 
oscillator, which we calculated to have an intrinsic hyper- 
polarizability of about 0.57. [lij Figure [5] shows the opti- 
mized non-perturbative hyperpolarizability when using a 
clipped harmonic oscillator as the starting potential. The 
properties of all of the optimized potentials are summa- 
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TABLE II: Summary of calculations with different starting 
potentials. f3 s is the hyperpolarizability of the starting po- 
tential while the other ones are after optimization. 



Function 


O " J \-J O f-^ is jt 


Pn p 


r (80) 
00 


Ar (80) 


V(x) 






(xlO -5 ) 


(xlO" 4 ) 





0.7089 0.7089 


0.7089 


37.8 


5.33 


30 tanh(a;) 


0.67 0.7084 0.6918 


0.7083 


779 


11.8 


X 


0.66 0.7088 0.7072 


0.7088 


78.7 


8.79 


x 2 


0.57 0.7089 0.7085 


0.7088 


18.6 


703 




0.68 0.7087 0.7049 


0.7087 


190 


9.76 


x + sin(x) 


0.67 0.7088 0.7073 


0.7088 


75.0 


8.46 


a; + losing) 0.04 0.7085 0.7085 0.7085 


1.65 


7.78 



> 





FIG. 5: Optimized potential energy function and first 15 
wavefunctions after 8,000 iterations. Starting potential is 
V(x) — x 2 , using the non-perturbative hyperpolarizability 
for optimization. 



rized in Table HU The clipped square root function also 
has a large hyperpolarizability (0.69). The optimized po- 
tential is shown in Figure^ In these cases, the amplitude 
of the wiggles are small and all the wavefunctions overlap. 
So, these fall into Class II. Note that the lack of wiggles 
shows that they are not an inevitable consequence of our 
numerical calculations. 

We may question whether small wiggles in the poten- 
tial energy function lead to large amplitude wiggles as an 
artifact of our numerical optimization technique. To test 
this hypothesis, we used the trial potential energy func- 
tion x + sin(x), where the wiggle amplitude is not large 
enough to cause the wavefunctions to localize at the min- 
ima. The optimized potential energy function retains an 
approximately linear from with only small fluctuation. 
In fact, the results are very similar to what we found for 
the linear starting potential and the wiggles do not affect 
the final result. The similarity between these cases can 
be seen in Table UU 

Next, we test a starting potential with large wiggles 
as shown in the upper portion of Figure [7J The lower 



FIG. 6: Optimized potential energy function and first 15 
wavefunctions after 8,000 iterations. Starting potential is 
V(x) = y^x, using the non-perturbative hyperpolarizability 
for optimization. 



energy eigenfunctions are found to be localized mostly in 
the first two wells. In fact, the lowest four energy eigen- 
functions are well approximated by harmonic oscillator 
wavefunctions, which are centrosymmctric. As a result, 
the first excited state holds most of the oscillator strength 
and the value of the intrinsic hyperpolarizability is only 
0.04. 

After 3000 interactions, this Class I potential energy 
function has high amplitude wiggles at a wavelength that 
is significantly shorter than the wavelength of the ini- 
tial sine function (bottom portion of Figure [7]). In com- 
mon with the optimized tanh(a;) function, the wiggles 
are of large but almost chaotically varying amplitude. 
This leads to wavefunctions that are spatially separated. 
While the wavefunctions are not as well separated as 
we find for the tanh(ai) starting potential, the optimized 
potential yields only two dominant transition from the 
ground state; so, this system is well approximated by 
a three-level model. As is apparent from Table HH the 
ground state sum rule (characterized by Tqq ^) is better 
obeyed in this optimized potential than in any others. 
So, the wavefunctions are accurate and all of the values 
of (3 have converged to the same value, suggesting that 
this calculation may be the most accurate of the set 

Our results bring up several interesting questions. 
First, all of our extensive numerical calculations, inde- 
pendent of the starting potential, yield an optimized in- 
trinsic hyperpolarizability with an upper bound of 0.71, 
which is about thirty percent lower than what the sum 
rules allow. Since numerical optimization can settle in to 
a local maximum, it is possible that all of the starting 
potentials are far from the global maximum of (3i n t = 1. 
Indeed, since most potentials lead to systems that require 
more than three dominant states to express the hyper- 
polarizability, this may in itself be an indicator that we 
are not at the fundamental limit precisely because these 
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FIG. 7: Potential energy function and first 15 wavefunctions 
before (top) and after (bottom) 3,000 iterations. Starting 
potential is of the form V(x) — x + 10sin(a;), using the non- 
perturbative hyperpolarizability for optimization. 

systems have more than three states. Indeed, the orig- 
inal results of Zhou and coworkers frames the problem 
in a way (i.e. a 15-level model in a potential limited to 
about 20 wiggles) that allows a solution to the optimiza- 
tion problem to lead to three dominant states. So, while 
it may be argued that this system is contrived and un- 
physical, we have found value in trying such toy models 
when testing various hypotheses. This toy model 

• leads to a three-level system as the three-level 
ansatz proposes 

• has the same qualitative properties as more precise 
methods 

• has given insights into making new molecules with 
record-breaking intrinsic hyperpolarizability 

Given the complexity of calculating nonlinear- 
susceptibilities, our semi-quantitative method may 
be a good way of generating new ideas. 



The three-level ansatz proposes that at the fundamen- 
tal limit, all transitions are negligible except between 
three dominant states. There appears to be no proof of 
the ansatz aside from the fact that it leads to an accurate 
prediction of the upper bound of nonlinear susceptibili- 
ties, both calculated and measured. To understand the 
motivation behind the ansatz, it is useful to understand 
how the two-level model optimizes the polarizability, a, 
without the need to rely on any assumptions. This is 
trivial to show by using the fact that the polarizability 
depends only on the positive-definite transition moments, 
(0| x \n) (n\ x |0), the same parameters that are found in 
the ground state sum rules. [2(| 

For nonlinear susceptibilities, the situation is much 
more complicated because the SOS expression depends 
on quantities such as (0| x \n) (n\ x \m) (m\ x |0), where 
these terms can be both positive and negative. Fur- 
thermore, the sum rules that relate excited states mo- 
ments to each other allow for these moments to be much 
larger than transition moments to the ground state. So, 
it would seem plausible that one could design a system 
with many excited states in a way that all of the tran- 
sition moments between excited states would add con- 
structively to yield a larger hyperpolarizability than what 
we calculate with the three-level ansatz. None of our 
numerical calculations, independent of the potential en- 
ergy function, yield a value greater than 0.71. Since our 
potential energy functions are general 1-dimensional po- 
tentials (i.e. the potentials are not limited to Coulomb 
potentials, nor are the wavefunctions approximated as is 
common in standard quantum chemical computations), 
our calculations most likely span a broader range of pos- 
sible wavefunctions leading to a larger variety of states 
that contribute to the hyperpolarizability. 

However, there appear to be local maxima associated 
with systems that behave as a three-level system and oth- 
ers with many states, and, the maximum values both are 
0.71. It is interesting that so may different sets of transi- 
tion moments and energies can yield the exact same local 
maximum. To gain a deeper appreciation of the under- 
lying physics, let's consider the transition moments and 
energies in the sum-over-states expression for the hyper- 
polarizability as adjustable parameters. For a system 
with N states, there are N — 1 energy parameters of the 
form E n — Eq. The moment matrix X{j has N 2 com- 
ponents. If the matrix is real, there are (iV 2 — N)/2 
unique off-diagonal terms and N diagonal dipole mo- 
ments. Since all dipole moments appear as differences 
of the form x nn — xoo , there are only N — 1 dipole mo- 
ment parameters. Therefore, the dipole matrix is char- 
acterized by (N 2 - N)/2 + N - 1 = (N + 2)(N - l)/2 
parameters. Combining the energy and dipole matrix pa- 
rameters, there are a total of (N + 2)(N - l)/2 + N - 1 
parameters. 

The N-state sum rules are of the form: 

(E n - X - {E m + E p )\ (m\ x \n) (n\ x \p) (7) 

71=0 ^ ' 



7 



h 2 N x 

~ 2m m *' 

so the sum rules comprise a total of JV 2 equations (i.e. an 
equation for each (m,p)). If the sum rules are truncated 
to JV states, the sum rule indexed by (m = JV, p = JV) 
is nonsensical because it contradicts the other sum rules. 
Furthermore, if the transition moments are real, then 
x mp = x pm , so only (JV 2 — N)/2 of the equations are 
independent. As such, there are a total of (JV 2 — N)/2 + 
N — 1 = (N + 2)(N — l)/2 independent equations. 

Since the SOS expression for the nonlinear- 
susceptibility has (JV + 2)(JV — l)/2 + JV — 1 parameters 
and the sum rules provide (JV + 2)(JV — l)/2 equations, 
the SOS expression can be reduced to a form with JV — 1 
parameters. For example, the three-level model for the 
hyperpolarizability, which is expressed in terms of 7 
parameters, can be reduced to two parameters using 5 
sum rule equations. In practice, however, even fewer 
sum rule equations are usually available because some 
of them lead to physically unreasonable consequences. 
While the (JV, JV) sum rule is clearly unphysical due 
to truncation, sum rule equations that are near equa- 
tion (JV, JV) may also be unphysical. In the case of 
the three-level model, it is found that the equation 
(2, 1) allows for an infinite hyperpolarizability, so that 
equation is ignored on the grounds that it violates the 
principle of physical soundness. [l2l l25l. [2^| This leads to 
a hyperpolarizability in terms of 3 variables, which are 
chosen to be E w , E — E 10 /E 2 q, and X — xi /x^ Ax . 
The expression is then maximized with respect to the 
two parameters E and X, leaving the final result a 
function of E±o. 

We conclude that the SOS expression for the hyperpo- 
larizability can be expressed in terms of at least JV — 1 
parameters; so, it would appear that as more levels are 
included in the SOS expression, there are more free pa- 
rameters that can be varied without violating the sum 
rules. As JV — > oo, there are an infinite number of ad- 
justable parameters. So, it is indeed puzzling that the 
three-level ansatz yields a fundamental limit that is con- 
sistent with all of our calculations for a wide range of 
potentials, many of which have many excited states. It 
may be that we are only considering a small subset of 
potential energy functions; or, perhaps the expression 
for the hyperpolarizability depends on the parameters in 
such a way that large matrix elements contribute to the 
hyperpolarizability with alternating signs so that the big 
terms cancel. This is a puzzle that needs to be solved if 
we are to understand what makes f3 large. 

To investigate whether the limiting behavior is 
due to our use of 1-dimcnsional potentials, we have 
also optimized the intrinsic hyperpolarizability in two- 
dimensions. In this case, we focus on the largest ten- 
sor component, /3 XXX and describe the potential as a 
superposition of point charges. As described in the 
literature. [l9j we solve the two-dimensional Schrodinger 



eigenvalue problem, 

h 2 

V 2 * + = E^ , (8) 

2m 

for the lowest ten to 25 energy eigenstates, depending on 
the degree of convergence of the resulting intrinsic hyper- 
polarizability. We use the two-dimensional logarithmic 
Coulomb potential, which for k nuclei with charges q\e, 
. . . , qke located at points . . . , s^ k \ is given by 

V(s) = ^5>log|| S - S «||, (9) 

where L is a characteristic length. With L = 2A, the 
force due to a charge at distance 2A is the same as it 
would be for a 3D Coulomb potential. 

We discretize the eigenvalue problem given by Equa- 
tion [8] using a quadratic finite element method (2ll. |27| 
and solve the resulting matrix eigenvalue problem for the 
ten to 25 smallest energy eigenvalues and corresponding 
eigenvectors by the implicitly-restarted Arnoldi method 
[22| as implemented in ARPACK [28|]. Each eigenvector 
yields a wave function \I/ n corresponding to energy level 
E n . The moments 

/oo /*oo 
/ si* m (si,s 2 )* n (si,s 2 )dsids2 
-ocJ — oo 

are computed, and these and the energy levels E n are 
used to compute ft 

Figure [5] shows the intrinsic hyperpolarizability of a 
two-nucleus molecule plotted as a function of the distance 
between the two nuclei and nuclear charge q\. The total 
nuclear charge is q\ + q 2 = +e, and is expressed in units 
of the proton charge, e. Three extrema are observed. 
The positive peak parameters are Pint = 0.649 for q\ = 
0.58 and d — 4.36A. The negative one yields (3i nt = 
-0.649 for qi = 0.42 and d = 4.361. The local negative 
peak that extends past the graph on the right reaches its 
maximum magnitude of /3j„t = —0.405 at q\ = 2.959 and 
d = 2.QA. 

Applying numerical optimization to the intrinsic hy- 
perpolarizability using the charges and separation be- 
tween the nuclei as parameters, we get = 0.654 at 
d = 4.5394, and q\ — 0.430 when the starting param- 
eters are near the positive peak; and j3i n t = —0.651, 
d = 4.443A, and q\ — 0.572 when optimization gives 
the negative peak. The peak parameters are the same 
within roundoff errors when optimization or plotting is 
used, confirming that the optimization procedure yields 
the correct local extrema. 

Figure [9] shows the intrinsic hyperpolarizability of an 
octupolar-likc molecule made of three evenly-spaced nu- 
clei on a circle plotted as a function of the circle's diam- 
eter and charge fraction e (q = ee). The charge on each 
of the other nuclei is e(l — e)/2. The positive peak at 
e = 0.333 and diameter D = 6.9A has a hyperpolariz- 
ability [3i n t — 0.326, while f3i nt — —0.605 for a charge 
fraction e = 0.44 and a diameter D = 6.8A 
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FIG. 8: The intrinsic hyperpolarizability of two nuclei as a 
function of the distance between them and the charge of one 
nucleus, qi where qi + qi — +e. 
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FIG. 9: The intrinsic hyperpolarizability of three evenly- 
spaced nuclei on a circle as a function of the circle's diameter 
and the charge, e (in units of e), on one of the nuclei. The 
charge on each of the other nuclei is e(l — e)/2. 

When the positions and magnitudes of the three 
charges are allowed to move freely, the best intrinsic 
hyperpolarizability obtained using numerical optimiza- 
tion is /3i nt = 0.685 for charges located at f\ = (0,0). 
r 2 = (-4.871,0.331), and f 3 = (-9.571,-0.161); with 
charges qi = 0.43e, q 2 = 0.217e, and q 3 = 0.351e. There 
are only small differences in the optimized values of fli nt 
depending on the starting positions and charges; and the 
best results are for a "molecule" that is nearly linear 
along the x-direction. This is not surprising given that 
the rax-component of flint is the optimized quantity. 

The two-dimensional analysis illustrates that numer- 



ical optimization correctly identifies the local maxima 
(peaks and valleys) and that the magnitude of maximum 
intrinsic hyperpolarizability (0.65 vs 0.68) is close to the 
maximum we get for the one-dimensional optimization 
of the potential energy function (0.71). All computa- 
tions we have tried, including varying the potential en- 
ergy function in one dimension or moving around point 
charges in a plane all yield an intrinsic hyperpolarizabil- 
ity that is less than 0.71. 

An open question is the origin of the factor-of-thirty 
gap between the best molecules and the fundamental 
limit, which had remained firm for decades through the 
year 2006. Several of the common proposed explana- 
tions, such as vibronic dilution, have been eliminated. fl4j 
Perhaps it is not possible to make large-enough varia- 
tions of the potential energy function without making 
the molecule unstable. Or, perhaps there are subtle is- 
sues with electron correlation, which prevents electrons 
from responding to light with their full potential. The 
fact that the idea of modulation of conjugation has lead 
to a 50% increase over the long-standing ceiling - reduc- 
ing the gap to a factor of twenty - makes it a promising 
approach for potential further improvements. Continued 
theoretical scrutiny, coupled with experiment, will be re- 
quired to confirm the validity of our approach. 



IV. CONCLUSIONS 

There appear to be many potential energy functions 
that lead to an intrinsic hyperpolarizability that is near 
the fundamental limit. These separate into two broad 
classes: one in which wiggles in the potential energy 
function forces the eigenfunctions to be spatially sepa- 
rated and a second class of monotonically varying wave- 
functions with small or no wiggles that allow for many 
strongly overlapping wavefunctions. Interestingly, all 
these one-dimensional "molecules" have the same max- 
imal intrinsic hyperpolarizability of 0.71. It is puzzling 
that the three-level ansatz correctly predicts the funda- 
mental limit even when the ansatz does not apply. A 
second open question pertains to the origin of the long- 
standing factor of 30 gap between the fundamental limit 
and the best molecules. The idea of conjugation modu- 
lation may be one promising approach for making wiggly 
potential energy profiles that lead to molecules that fall 
into the gap. Given that there are so many choices of 
potential energy functions that lead to maximal intrinsic 
hyperpolarizability, it may be possible to engineer many 
new classes of exotic molecules with record intrinsic hy- 
perpolarizability. 
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